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Abstract. We describe a “spatio-spectral” deconvolution algorithm for 
wide-band imaging in radio interferometry. In contrast with the existing 
multi-frequency reconstruction algorithms, the proposed method does 
not rely on a model of the sky-brightness spectral distribution. This non- 
parametric approach can be of particular interest for the new generation 
of low frequency radiotelescopes. The proposed solution formalizes the 
reconstruction problem as a convex optimization problem with spatial 
and spectral regularizations. The efficiency of this approach has been 
already proven for narrow-band image reconstruction and the present 
contribution can be considered as its extension to the multi-frequency 
case. Because the number of frequency bands multiplies the size of the 
inverse problem, particular attention is devoted to the derivation of an 
iterative large scale optimization algorithm. It is shown that the main 
computational bottleneck of the approach, which lies in the resolution 
of a linear system, can be efficiently overcome by a fully parallel imple¬ 
mentation w.r.t. the frequencies, where each processor reconstructs a 
narrow-band image. All the other optimization steps are extremely fast. 

A parallel implementation of the algorithm in Julia is publicly available 
at https : //github. com/andf errari. Preliminary simulations illustrate 
the performances of the method and its ability to reconstruct complex 
spatio-spectral structures. 
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1. Introduction 


Recently, much attention has been paid to the development of image re¬ 
construction algorithms for the incoming and future radio facilities. Most 


recent contributions, such as ( 

Carril 

o et al. 

2012 

Dabbech et al. 

2015 

; Car- 

rillo et al. 

2014 

) or (Garsden et al., 

2015) heavily rely on sparse estimation. 

In the wake of CLEAN algorithm ( 

Hogbom, 1974) and its multiresolution 


variants, sparse models have indeed proved in the last decades to be a pow- 
erfull approach for radiointerferometric image reconstruction in particular, 
and for the resolution of inverse problems in general. 

In addition to their high spatial resolution, the wide bandwidths of the 
new generation of radio interferometers makes possible the reconstruction 
of complex spectral structures. The Square Kilometre Array (SKA) and 
its precursors will achieve (sub-)arcsec resolution over hundreds of MHz in¬ 
stantaneous bandwidths and with a tremendously broad band coverage (see 


Table 1 in the SKA1 System Baseline Design document (Dewdney et al. 


2013)). The reconstruction of both spatial and spectral behaviour of con¬ 


tinuum radio sources is an essential tool to characterize the astrophysical 


origin of their detected radiation, e.g. (Kraus, 1986). (Rau and Cornwell 


2011) opened the way to multi-frequency deconvolution algorithms, which 


aim to reconstruct simultaneously spatial and spectral structures. The ap¬ 


proach proposed in (Rau and Cornwell, 2011) relies on the parameterization 


of the frequency-dependent brightness distribution as a power law with a 
varying index. A Taylor expansion is adopted to model the flux dependence 
in frequency of astrophysical radio sources, whose synchrotron or thermal 
spectra can be described by power-laws (Conway et al., 1990). From an 


estimation point of view, the ratio of the additional unknowns introduced 
by such a multi-frequency model (e.g. the average spectral-indexes and the 
spectral-curvatures for a second order model) to the additional equations 
(the multiplication of measurements by the number of frequency bands) is 
clearly in favor of a multi-frequency reconstruction approach. Explicit mod¬ 
els of the frequency dependence of radio sources have also been introduced 


and exploited in (Junklewitz et al., 2014) and (Bajkova and Pushkarev 


2011). In (Junklewitz et al., 2014), the authors propose to address the es¬ 


timation problem using a Bayesian framework. (Bajkova and Pushkarev 


2011) proposes to constrain a maximum entropy estimation algorithm in 


order to account for the frequency dependence of the intensities. 

These “semi-parametric” methods rely on spectral models and thus clearly 
offer advantages and estimation accuracy when the model is indeed appro¬ 
priate. Across the broad frequency coverage of current radio facilities how¬ 
ever, radio sources exhibiting complex spectral shapes (not simple power 
laws) can be expected. Such sources can show one or more relative min¬ 
ima, breaks and turnovers ( |Kellermann , 1974). More recently, (Scaife and 
|Heald , |2012 ) have shown that second order broadband spectral models are 
often insufficient for the new generation of low frequency telescopes such 
as the Low Frequency Array (LOFAR). A non-parametric approach in the 
multi-frequency reconstruction of radio sources is also definitely needed for 
the full Stokes (Q, U, V) wide-band imaging, where Taylor expansion is not 
the appropriate physical model to be adopted (Rau and Cornwell, 2011). In 
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the same direction, (Wenger and Magnor, 2014) have proposed to relax the 


spectral power-law model by formulating the problem as an inverse problem 
with a smooth spectral regularization allowing for local deviations. 

This communication proposes to fully formalize multi-frequency image 
reconstruction as a regularized inverse problem. As such, it extends the 


models proposed in (Carrillo et al., 2012, 2014; Garsden et al., 2015) to the 


multi-frequency observation mode by adding a spectral regularization term. 
It is worthy to note that the proposed approach shares formal similarities 
with a spatio-spectral image reconstruction algorithm recently proposed for 


optical interferometry (Schutz et al., 2014; Thiebaut et al., 2013). 


Section 2. introduces the data model and the inverse problem criterion. 
Section 3. derives an efficient optimization algorithm to minimize the re¬ 
lated convex problem. The time consuming steps of this iterative algorithm 
can be parallelized frequency-wise leading to an overall computation time 
comparable to a narrow-band reconstruction algorithm. Section 3. presents 
preliminary simulation results where the observations are obtained using the 


MeqTrees package (Noordam and Smirnov, 2010). 


2. Models and notations 

We denote by X£ the column vector collecting the sky intensity image at 
observing frequency v#, £ = 1 ,... ,L. The sky image map X£ is related to 
the “dirty image” at frequency z/^, denoted as by 


y e = H e x e + ri£ 


( 1 ) 


where ri£ is a noise vector. Hi is a convolution matrix containing the single¬ 
frequency point-spread function (PSF), which includes weighting factors and 


visibility mapping on a spatial frequency grid (Rau and Cornwell, 2011). 


Stacking the images at different frequencies in a single vector leads to the 
wideband model 


y Hx + n 


( 2 ) 


where H is a block diagonal matrix containing the Hi matrices on the main 
diagonal. 

Note that we consider herein without loss of generality the model con¬ 
necting the sky intensity to the dirty map. The same algorithm as the one 
described below can be derived from the model connecting the sky intensity 
to the complex visibilities by replacing in ([!]) yi by the sample visibilities 
and Hi by the operator mapping the sky intensity to the visibilities at the 
observing frequency z/^, see Appendix A.2 of (Rau and Cornwell, 20lT|. 

Eq. [2] defines an ill-posed linear inverse problem. Among all existing 
methods to solve inverse problems, we will focus on sparse regularization 
and optimization techniques, which have gained much popularity during 
the last decade. The regularization game consists in minimizing a criterion 
composed by a fidelity term and a regularization term / reg linked to some 
prior on the solution. We will consider here the objective: 

min ^\\y - Hx\\ 2 + f reg (x) 
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Many recent works have shown that regularization based on sparse repre¬ 
sentations in appropriate transform domain(s) of the intensity map being re¬ 
stored can be very effective. This approach can be formulated in an analysis 
framework: the criterion is minimized directly w.r.t. the image parameters, 
i.e., x. An alternative formulation, not considered in the present study, is 
the synthesis framework: the criterion is then minimized w.r.t. the coeffi¬ 
cients of the image decomposition. Those can be much more numerous than 
the image parameters if the decomposition is redundant - and it should to be 
efficient. This leads to an increased computational load w.r.t. the analysis 
framework. These two formalisms are discussed in (Elad et ah, 2007). In 
practice and performance-wise, it is still unclear which approach should be 
preferred and under which conditions. The efficiency of both approaches for 


narrow-band radio-interferometric imaging has been proven in Carrillo et al. 
(|2012[ 2014) (analysis), Dabbech et al. (2015) (analysis-synthesis), |Garsden 


et al. (2015) (synthesis). Similarly to the first two references, we will focus 


on a sparse analysis prior and in complement to the classical positivity and 
quadratic smoothing regularization, we will consider a regularization of the 
form: 


f re %(x) = 


-i . ( 


Ur . 


where X is the matrix defined as X = (aq, ..., xi u ) and W s and W v are the 
matrices associated with respectively the spatial and spectral analysis regu¬ 
larizations. In practice these matrices take the form of (usually redundant) 
dictionaries with dedicated fast transforms, some examples will be given be¬ 
low. It is important to underline the central role of the last regularization 
term with parameter in Eq. [4j This term prevents the optimization prob¬ 
lem <© to be separable w.r.t. the X£. This makes the sparse spatial and 
spectral priors imbedded and the regularization truly spatio-spectral. 


3. Optimization algorithm 


We propose to find the solution of the convex problem (3|4) using the 
Alternating Direction Method of Multipliers (ADMM) algorithm. For a 


recent comprehensive review of ADMM see (Boydetahl 2011). Convergence 


of ADMM was demonstrated in ( |Eckstein and Bertsekas , 1992). As it will 
be shown below, this method is particularly interesting to solve large-scale 
problems such as as it leads to successive steps that can be parallelized 

w.r.t. the images at each frequency (i.e., the aq, columns of X) or w.r.t. 
the spectra at each pixel position (i.e., the rows of X). 

The end of this section is devoted to the derivation of this algorithm. The 
optimization problem ( 3j4 ) is equivalent to: 


min ^\\Y-HX\\ 2 + 1 r+ (P) + ^\\X\\ 2 f + i, s \\T\\ 1 +i, u \\V\\ 1 (5) 
subject to: P = X, T = W S X , V = SW U , S = X (6) 


where Y = (j/ 1? ..., y^). Auxiliary variables P, T, V and S have associ¬ 
ated Lagrange multipliers r p , T t , T v and V s and augmented Lagrangian 
parameters pp, pp, py and ps- Denoting with a + superscript the updated 
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quantities, the alternated minimizations of the augmented Lagrangian of 


a 


give: 


I. Minimization w.r.t. X. This step operates separately on each fre¬ 
quency image X£. It requires to solve for each frequency t/£ the linear 
system 


^ - 1 . . . , Ijy 


( 7 ) 


where: 


Q e = HjH e + p T WjW s + (^ + p s + pp)I (8) 

b = h t y + wj ( r T + Pt t) + r p + Pp p + t s + Ps s ( 9 ) 

and b£ are the columns of B. 

II. Minimization w.r.t. P. Defining P — X — Pp 1 Tp, this minimization 
simplifies to the proximity operator: 

min 1 R+ (P) + ^||P-P||^ (10) 

and leads to the positive projection of each element Pij: 

Pfj = m ax(0, Pij) (11) 

III. Minimization w.r.t. T. Defining T = W S X — minimisation 

w.r.t. T simplifies to the proximity operator: 

min ^Ill’ll! + ^||T-T||2. (12) 

T A 

Consequently, each element T^j of T is updated according to the 
soft thresholding operator: 

Tt d =Tij ma^|o,l-|^J (13) 

IV. Minimization w.r.t. V. Defining V = SW U — p^Ty, minimisation 
w.r.t. V simplifies to: 

mm Pl/ \\V\\ 1 + ^\\V-V\\ 2 F (14) 

Each element Vij of V is updated according to the soft thresh¬ 
olding operator: 

V tj = Vij max 1 - ( 15 ) 

V. Minimization w.r.t. S. This steps operates separately on each voxel. 
It requires to solve the linear system S+R = C where: 

R = P yW v Wl + p s I (16) 

C = (T V + Pv V)Wj -T s + PS X (17) 

This step is the twin of step I. (minimization w.r.t. X). Choosing 
for W v an orthonormal transform, or a union of m v orthonormal 
transforms, is of particular interest as it leads to: 

- 1 —— ((T v + Pv V)W T v -T s + PS X) 

m uPv + V / 
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VI. Update of the Lagrangian multipliers. Finally, the Lagrangian mul¬ 
tipliers are updated according to the standard way: 


r+ = r P + Pp (p - x) 

(19) 

r+ = r T + Pt (t - w s x ) 

(20) 

T+ = T v + p v (V-SW u ) 

(21) 

r£ = r s + PS (s - x) 

(22) 


The previous six steps are iterated until convergence (see (Boyd et al. ,2011) 
for the stopping criterion). This leads to algorithm [l] where S r (-) denotes 
the element-wise soft thresholding operator with threshold r. 


Initialize to zero X , P, T, V , S', Tp, Tp, Ty and Tp; 

repeat 

/* update primal variables */ 

for t — 1..., l v do in parallel 

Compute b£ using Q 

/* conjugate gradient algorithm */ 

Solve = bi where is given by Q 

end 

P max(0, X — pp 1 Tp) 

T^S^ s/pt (W s X-pPt t ) 

V^S^ /pv {SW„- P y l T v ) 

s <- (m u pv + ps)~ l ((IV + PvV)Wl - Vs + p S X) 

/* update dual variables */ 

T P ^Tp + p P (P - X ) 

r T «- r T + Pt {t - w s x) 

T v ^T v + pv(V -SWu) 

Ts^T s + Ps(S-X) 
until stopping criterion is satisfied. 

return X 

Algorithm 1: Multi-frequency reconstruction alorithm. 


Similarly to the narrowband case (Carrillo et al., 2014), resolution of 0 
in step I. is the bottleneck of the algorithm. Note that the possibility to 
solve the l v systems in parallel does not increase the computation time of 
this step compared to the narrowband case. 

As in step V., computation of the second term of Qf defined in Eq. H 
simplifies when W s is a union of n s orthonormal bases: 

Qi — Hj Hi + pi, p = p e + Ps + Pp + nspT (23) 

The resolution of each linear system can then be drastically accelerated using 


a conjugate gradient algorithm (Hestenes and Stiefel, 1952). However, it is 


worth noting that in the dirty image model ([]J) Hi is a convolution matrix 
and consequently Qf in (23) is also a convolution matrix. Therefore Qf X bi 
can be efficiently computed by the convolution of bf with the filter with 
frequency response (\hi(m, n)\ 2 + p) _1 where hi is the Fourier transform of 
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the point spread function associated to H£. Note the obvious role of p as a 
regularizer of this inversion when expressed in the frequency domain. 


4. Simulations 


A Julia implementation of the code is availabl<0 It takes the PSFs and 
dirty images for entries. W s is a dictionary composed by the concatena¬ 
tion of the first eight orthonormal Daubechies wavelet bases (Dbl-Db8) as in 
(Carrillo et ah, 2014) and a Haar wavelet basis. In order to promote smooth¬ 


ness W v implements a Discrete Cosine Transform (DCT). As in (Wenger 


and Magnor, 2014), Chebyshev polynomial basis functions should be also 


appropriate. The implementation takes advantage of the Julia language 
parallel processing support. 

The PSFs used in this simulation where computed using the MeqTrees 
package (Noordam and Smirnov, 2010) with MeerKAT arrays configuration. 
We simulated t v — 15 frequency bands. The first is centered at 1.025 GHz. 
The central frequencies of the 15 bands are separated by 50 MHz. The total 
observation time is 8 hours. The size in pixels of the images are 256 x 256. 
Figure [l] shows 50 x 50 pixels images of the center of the PSF at four fre¬ 
quencies (1.025 HZ, 1.225 GHz, 1.475 GHz and 1.725GHz). The parameters 
p e , p s and fi v control the levels of regularization (resp. smoothing, spatial 
and spectral). They have been set to p e = 0.001, p s = 0.5 and ji v — 1.0. 
The augmented Lagragian parameters pp, pp, py and ps affect the speed 
of convergence. They have been set to pp = 1.0, pp = 5.0, py = 2.0 and 
ps = 1.0. The classical termination criterion for the ADMM algorithm is 
that the primal and dual residuals must be small (Boyd et ah, 2011). In 


order to avoid additional parameters, the algorithm will be simply stopped 
after 350 iterations in the following simulations. 

The first simulation evaluates the ability of the algorithm to reconstruct 


spectra. The sky image is similar to the first simulation of (Rau and Corn¬ 


well, 2011). It consists in two overlapping Gaussian profiles centered at pixel 


108 and 148, with constant spectral indices equal to respectively -1.0 and 


+ 1 . 0 . 

Figure [2] shows in four frequency bands: the sky image, the dirty im¬ 
age and the reconstructed sky obtained using algorithm [lj Figure [3] shows 
the spectra associated to the pixel at the center of each Gaussian. In each 
case the theoretical spectrum and the reconstructed spectrum are plotted. 
In order to derive a synchrotron spectral index map, each spectrum of the 
original sky (which equals here a linear combination of v and 1/z/) is fitted 
to a second order power law model; i.e. oTog(z/) + (3 log 2 (z/) in a log scale. 
Spectral indices are then estimated fitting the same model to the recon¬ 
structed spectra. Figure [4] shows the estimated a for the pixels on the line 
joining the center of the two Gaussian profiles. The spectral indices esti¬ 
mated directly from the dirty image are given in the figure for comparison. 
These preliminary simulations show the ability of the proposed method to 
reconstruct complex spectral signatures using a non parametric approach. 

The scope of the second simulation is to evaluate the performances of the 
method in a more realistic scenario. Radio emission from an HII region in the 


1 https://github.com/andferrari/muffin.jl 
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v= 1.025 GHz v= 1.225 GHz 



Figure 1. PSFs of MeerKAT in 4 among the 15 simulated 
frequency bands. The PSFs are obtained from MeqTrees. 
The total observation time is 8 hours. The images are a 
50 x 50 pixels zoom of the central region. 


M31 galaxy is used as a reference sky image. A sky cube is then computed 
applying a first order power-law spectrum model to the M31 image. The 
256 x 256 map of spectral indices is constructed following the procedure 


detailed in (Junklewitz et ah, 2014): for each pixel the spectral index a is 


a linear combination of an homogeneous Gaussian field and the reference 
sky image. The first correlates the spectral indices spatially and the second 
correlates in a pixel the spectral index with the object intensity. Figure [5] 
shows the sky map and a typical spectral indices map. A white gaussian 
noise with a constant variance corresponding to a signal to noise ratio of 20 
dB has been added to all the dirty images. 

The computation time required to reconstruct the 256 x 256 x 15 pixels 
cube, using — 15 core^] is approximately 1.5 hours. Figure [6] shows in 
four frequency bands: the sky image, the dirty image and the reconstructed 
map. Figure [ 7 ] shows for the same frequency bands the error images between 
the sky and the reconstructed sky (square root of the absolute value of the 
difference between the two images). Finally, figure [8] gives the relative root 
mean square error (RMSE) as a function of the frequency band. These 
results clearly show the capacity of the algorithm to recover details of the 
sky image at every frequency. 


5. Conclusion 

This paper presents first results on multi-frequency image reconstruction 
for radio interferometry using a fully regularized inverse problem approach. 
The proposed algorithm relies on the minimisation of a data fidelity term 

2 http://calculs.unice.fr 
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Figure 2. Multi-frequency reconstruction results. Each line 
is associated to a frequency band. The first column is the 
original object. The second column is the dirty image and the 
third column the reconstructed object. The second column 
(dirty images) is scaled down in order to fit to the color map 
of the two other columns. 

spatially and spectrally regularized. An efficient iterative algorithm is de¬ 
rived to minimize the related convex cost function. The computational 
bottleneck of the algorithm can be parallelized w.r.t. the images at each 
frequency leading to an overall computation time of the order of narrow- 
band algorithms. Future work will focus among others on the derivation of 
specific spectral regularizations and on the necessity of taking into account 
frequency-dependent instrumental effects such as a the primary-beam. 
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Figure 3. Spectra associated to the pixel at the center of 
each Gaussian profile. Sky model in blue and reconstructed 
sky model in green. 



pixel 


Figure 4. Estimated spectral indices for the pixels on the 
line joining the center of the two Gaussian profiles. Sky 
model in blue, reconstructed sky model in green and dirty 
image in red. 
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Figure 5. Radio emission from an HII region in the M31 
galaxy and a spectral indices map generated combining the 
sky intensity and a Gaussian homogeneous random field. 
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Figure 6. Multi-frequency reconstruction results. Each 
line is associated to a frequency band. The first column is 
the original sky. The second column is the dirty image and 
the third column the reconstructed sky. The second column 
(dirty images) is scaled down in order to fit to the color map 
of the two other columns. 
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Figure 7. Error between the sky and the reconstructed 
sky. The images show the square root of the absolute value 
of the difference between the two images. 



Figure 8. Relative root mean square error (RMSE) of the 
reconstructed sky as a function of the frequency. 
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